function [RFE,convKernels]=timeInvariantMultiCell(lfp,spikes,KernelWidth,penalty)
[N,numCells]=size(spikes);
if nargin<4
    penalty=1e-16;
elseif penalty<1e-16
    penalty=1e-16;
end
%% Initialization
% Get initial values
for n=1:numCells
    [convKernels(:,n)]=getBoundedWiener(lfp,spikes(:,n),KernelWidth,penalty);
end
% Make predictions from inital value
lfpPred=zeros(numel(lfp),numCells);
for n=1:numCells
    tmp=conv(spikes(:,n),convKernels(:,n),'full');
    lfpPred(:,n)=tmp(KernelWidth+1:end-KernelWidth);
end
a=lfpPred\lfp;
lfpPred=bsxfun(@times,a(:)',lfpPred);
lfpHat=sum(lfpPred,2);
convKernels=bsxfun(@times,a(:)',convKernels);
RFE=1-sum((lfp-lfpHat).^2)./sum(lfp.^2);
fprintf('RFE after initialization: %0.4f\n',RFE);
%% Fine-tuning:
for iter=1:5
    for n=1:numCells
        lfpHat=lfpHat-lfpPred(:,n);
        [convKernels(:,n)]=getBoundedWiener(lfp-lfpHat,spikes(:,n),...
            KernelWidth,penalty);
        tmp=conv(spikes(:,n),convKernels(:,n),'full');
        lfpPred(:,n)=tmp(KernelWidth+1:end-KernelWidth);
        lfpHat=lfpHat+lfpPred(:,n);
    end
    RFE=1-sum((lfp-lfpHat).^2)./sum(lfp.^2);
    fprintf('RFE after %d fine-tunes: %0.4f\n',iter,RFE);
end
